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ABSTRACT 

We study in detail the motions of three planets interacting with each other 
under the influence of a central star. It is known that the system with more than 
two planets becomes unstable after remaining quasi-stable for long times, leading 
to highly eccentric orbital motions or ejections of some of the planets. In this 
paper, we are concerned with the underlying physics for this quasi-stability as 
well as the subsequent instability and advocate the so-called "stagnant motion" 
in the phase space, which has been explored in the field of dynamical system. We 
employ the Lyapunov exponent, the power spectra of orbital elements and the 
distribution of the durations of quasi-stable motions to analyze the phase space 
structure of the three-planet system, the simplest and hopefully representative 
one that shows the instability. We find from the Lyapunov exponent that the 
system is almost non-chaotic in the initial quasi-stable state whereas it becomes 
intermittently chaotic thereafter. The non-chaotic motions produce the hori- 
zontal dense band in the action-angle plot whereas the voids correspond to the 
chaotic motions. We obtain power laws for the power spectra of orbital eccen- 
tricities. Power-law distributions are also found for the durations of quasi-stable 
states. All these results combined together, we may reach the following picture: 
the phase space consists of the so-called KAM tori surrounded by satellite tori 
and imbedded in the chaotic sea. The satellite tori have a self-similar distribu- 
tion and are responsible for the scale-free power-law distributions of the duration 
times. The system is trapped around one of the KAM torus and the satellites for 
a long time (the stagnant motion) and moves to another KAM torus with its own 
satellites from time to time, corresponding to the intermittent chaotic behaviors. 

Subject headings: celestial mechanics — planetary systems — solar system: general 
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Introduction 



■I 



More than 300 exoplanets have been discovered so far |j and, interestingly, some of them 
have a quite different appearance from that of our solar system. The existence of so-called 
"eccentric planets", that is, pla nets with high orb ital eccentricities, for example, attracts 



attentions of many researchers ( iButler et al. 



20061 ). Stimulated by these observations, the 



planetary formation theo ry including the orig in of the eccentric planets has made substantial 



progress over the years (iKokubo et al. 



20061 ). In the standard theory, terrestrial planets 



are thought to be formed through giant impacts of prq to-planets (planets with sub-Earth 



masses) in crossing orbits ( IChambers fc Wetherill 



19981 ). Since N-body simulations suggest 



that proto-planets ar e formed in nearly circular orbits separated by several Hill radii 



( iKokubo &: Ida 



19951 ). some destabilizing processes are expected to operate and make 
proto-planets originally in the well-separated circular orbits collide with each other and 
grow up to terrestrial planets. 

The stability of the system with two planets around a central star has been thoroughly 
investigated in celestial mechanics and it is known that there exists a critical orbital 
separation between the planets, beyond which the planets never experience close encounters 



and the system remains stable forever (IMarchal fc Bozis 



1980 



Gladman 



19931). The 



situation changes dras tically, however, i 



numerical simulations. 



Chambers et al. 



anot her planet is added to the system. Using 



(119961 ) demonstrated that the systems with more 



than two planets become unstable even for large orbital separations. Although the planetary 
motions remain regular for some time at first, one of the planets eventually comes close 
enough to another (that is, within the Hill radius of the latter), leading to subsequent 
orbital crossings. 



^http:/ /exoplanet.eu for the latest information. 
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This is a good news for the terrestrial 
formation of the eccentric planets. In fact, 



ormation theory and may a. 



Marzari fc Weidenschilling ( 12002| ) numerically 



so account for the 



integrated the motions of three Jupiter-mass planets and found in most of their simulations 
that one of the planets is ejected from the system and the others are left in the sys tem with 



high eccentricities. Several attempts ( lJuric fc Tremaine 



2007 



Ford fc Rasio 



20071 ) have 



been made to reproduce the observed eccentricity distribution by the orbital instability. 
Their results seem to be consistent wit h the observations although the latter itself may be 



somewhat biased (jShen fc Turner 



20081 ) 



The planetary motions in these systems are interesting in their own right. As 
mentioned above, we commonly observe a long period of quasi-regular motions that look 
like independent Keplerian motions before the eventual orbital crossings. The switch is 
sudden and quick. The duration of t he quasi- regular motions is sensitive to the i nitial 



orbital separations and eccentricities (j Chambers et al 



1996 



Yoshinaga et al 



Zhou et al 



19991). Using 



(120071 ) claimed that the 



numerical simulations and simplified analytical models, 
gradual deviation from the Keplerian motions can be regarded as a random walk process. 
These efforts notwithstanding, the underlying physics behind the phenomena such as the 
long period of quasi-regular motions followed by the sudden transition to chaotic states is 
remaining to be revealed and is the main concern of this paper. 

We attempt to understand this phenomenon as the so-called "stagnant motion" in the 
phase space, which will be described below. We pay attention to a similar phenomenon 
known in the field of dynamical system. In nearly integrable systems, a sudden transition 
from a regular motion sustained for a long time to chaotic motions is often observed. In 



the nonlinear lattice problems, for example. 



Hirooka fc Saitd (119681 ) found that an initially 



imposed normal mode experiences sudden energy exchanges among several other modes 
after long regular oscillations. They called this "the induction phenomenon" and referred 



- 5 - 



to the duration of the regular motion as "the induction period" (jSaito et aLlll970l ) 



Aizawa et al 



( 119891 ) constructed a so-called "stagnant motion" model 



phen omenon. According to the KAM theorem (IKolmogorov 



1979 



Moser 



'or this induction 



1958 



Arnol'd 



19631 ). the phase space of a nearly integrable system retains tori, which exist in the integrable 
system, if perturbations to the integrable system are sufficiently weak. It is generally 
expected that the so-called KAM tori will survive even for not so small perturbations. In 
the stagnant motion model, it is assumed that the KAM torus exists in "the chaotic sea", 
the region corresponding to chaotic motions of the system, being surrounded by a thin layer 
called "the stagnant layer", in which smaller tori are distributed in a self-similar manner 
(see Fig. [1]). The system shows nearly regular behavior when the phase space orbit is 
trapped in the stagnant layer whereas the sudden transition to chaotic motions occur when 
the orbit escapes out of the layer. This model is successful in reproducing the scale-free 
power spectra and the distribution of the induction period. In this paper we show some 
evidence to support the interpretation of the motions of three-planet system as one of the 
induction phenomena and attempt to understand them in the frame work of the stagnant 
motion model. 



The organization of the paper is as follows. We summarize the numerical models in 
section 2. In section 3, we describe the methods of analysis. The results are presented in 
section 4 and the summary and discussions are given in section 5. 



2. Models 

In this paper, we restrict the investigation to the simplest multi-planet system, which 
shows the behavior mentioned above: the system consists of a central star with M^, = IMq 
and three planets with an identical mass in coplanar orbits. We consider two cases for the 
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planetary mass, rripi: (1) rripi = 10"''' Mq (the proto-planet system) and (2) rripi = 10~^Mq 
(the Jupiter system). The initial semi-major axis of the innermost planet is set to be 
ai = lAU for the prot o-planet system and ai = 5AU for the Jupiter system. Following 



Chambers et al. 



( 1l996l ). we give the initial radial locations of the outer planets as 

Qi+i = CLi + Ar/i(j^j+i), ( 

where Oj is the semi-major axis of the i-th planet counted from the innermost one, and 
i^h{i,i+i) is the Hill radius for the pair of i-th and i + 1-th planets defined as 

I \ 1/3 I 

rrii + mj\ Oj + aj 



-Mm)=(^^ ^— (2) 



In addition, we also consider the case with rripi = 10 ^Mq, ai = 5AU, a? = 7.25AU, 



g.-^ = 9.5AU, that is, the same parameter set as that used in 



Marzari fc Weidenschilling 



( I2OO2I ) except for no inclination in our model. 

The Hamiltonian in the barycentric coordinates consists of three parts and is given as 

H = Hxep + -f/^* + Hint, (3) 

1=1 ^ ' 

2mo / ' 

In the above equations, G is the gravitational constant and mj, rj denote the mass, 
momentum and radial position of the i-th object, where i = corresponds to the central 
star. Tij stands for the distance between the ^-th and j-th objects. Hxep is an integrable 
Hamiltonian corresponding to independent Keplerian motions of three planets with respect 
to the barycenter. H^, is a correction originated from the orbital motion of the central star 
itself by the attraction of the planets. Hint expresses the interactions between the planets. 
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For each model listed in Table [H we generate 10000 different initial conditions, which 
have the same integrals of motion, that is, the linear and angular momenta and total energy. 
This constraint is important to explore the structure in the phase space. The orbital phases 
of planets are given randomly. Since the radial location of the innermost planet, ai, is 
fixed, the remaining parameters, the radial locations of the outer planets, a2, 03, and the 
velocity of the central star, Vsx, Vgy, are determined so that the system should have the 
same values of the integrals of motion. As a result of this constraint, the resultant initial 
radial locations of the planets differ only slightly among 10000 realizations. The planets 
have slight eccentricity ( e ~ 10"'' for rripi = 10~^ and e ~ 10"^ for rripi = 10~^) initially 
because of the non-zero velocity of the central star relative to the barycenter. Using the 
ensemble obtained in this way, we obtain various distributions and take statistics thereof. 
We summarize in Table [T] the input parameters (ai and the integrals of motion) as well 
as the orbital separations and eccentricities averaged over the 10000 realizations for each 
model. 

Numeric al integrations are performed with the MERCURY6 package, which was 



developed by IChambersI ( I2OOOI ). For the proto-planet system, we integrate the orbital 
motions until the first close encounter occurs. For arbitrarily chosen three models among 
10000 realizations, we continue the integration after the close encounter up to lO'' years in 
order to compute the Lyapunov exponent and power spectra in the post-encounter phase. 
For the Jupiter system. The integration is terminated when one of the planets is ejected 
from the system. 

Before closing this section, we discuss the relative magnitude of each part of the 
Hamiltonian given in Eq. (|3]) and the existence of KAM tori in our system. The relative 
magnitude of if* to Hxep is always of the order of rripi/M^, and so is the ratio of Hint to 
Hxep unless Vij/viQ becomes as small as ~ mpi/M^. In our models, the value of mpi/M^ 



- 8 - 



is 10"'' for the proto-planet system or 10~^ for the Jupiter system and the minimum 
value of Tjj/rjo, which is achieved when two planets have the same orbital phase, is 
~ {fTT'pi/M^:Y^^ ^ nipi/M^:. Thus, and Hint are always small in our models. 

As mentioned in Introduction, the KAM tori exist in the phase space of the nearly 
integrable system, whose Hamiltonian is expressed as the sum of the integrable part Hq{I) 
and the perturbation eH{I, 9): 

H{I,e) = Ho{I) + eH{I,d) (e«l), (4) 

provided the perturbation is sufficiently small. Here / and 6 are the action and angle 
variable, respectively. It should be noted that our Hamiltonian does not meet this condition. 
Hence we look for evidence that this is really the case, employing the Lyapunov exponent, 
power spectra of orbital elements and induction periods, in addition to the trajectories in 
the 1 — 6 plane, which will be described in the following section. 

3. Analysis Methods 

In order to get some insight into the phase space structure of our models, we employ 
three tools, that is, the Lyapunov exponent, the power spectra of orbital elements and the 
distribution of the duration of the quasi-regular motions. The first two are useful to see 
the degree of chaos of the system. If there remain KAM tori in the phase space indeed, 
the system is expected to show both non-chaotic and chaotic features, which will then 
be captured by these measures. The last quantity will tell us if the planetary motions of 
our models can be interpreted as the induction phenomenon. In fact, the distribution is 
expected to have a power law if it is really the case. 
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3.1. Lyapunov exponent 

We compute the so-called maximum Lyapunov exponent. The maximum global 
Lyapunov exponent, Xgiobah is the local growth rate of the distance, || (5x(t) ||, between 
adjacent orbits in the phase space and is defined more precisely as 

, - r ^ II S^jT) II 

A,,„,.. = ^lim -logp^. (5) 

It is, of course, impossible in practice to compute the growth of the distance over the infinite 
time. We calculate instead the following quantity as a function of time and look into their 
behavior: 

1 II 5x(t) II 

A,.o..Kt) = -logp^. (6) 
The asymptotic limit of this function at t — )■ oo gives the original Lyapunov exponent. 

The integrable system has a null Lyapunov exponent while chaotic systems have a 
positive Lyapunov exponent. If the system is nearly integrable in particular, the Lyapunov 
exponent defined by Eq. ([6]) oscillates around a small but finite positive value and does not 
converge as t — )■ cxd. 

In analyzing the system that shows both quasi-regular and chaotic behaviors 
alternatively, it is useful to look also into the local Lyapunov exponent, Xiocah which is the 
same local growth rate of the orbital separation in much shorter times and is expressed as 

/ \ 1 , II ^xfrir) II 

Xiocaiin, T} = - log — TV^' in 

T II ox((n — l)r) II 

where n specifies an interval with a period of r. The choice of r is rather arbitrary. It 
should be longer than the typical orbital period but shorter than the time scale of the 
quasi-regular or chaotic motions of interest. If chosen appropriately, it will indicate the 
local degree of chaos. 

In the following, the interval, r, for the local Lyapunov exponent is chosen to be 100 
years, which corresponds to 10 ~ 100 times the orbital periods. Both Xgiobai and Xiocai are 
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obtained by numerically integrating the linearized equations of motion along th e phase 



space orbit given by the integration of the equations of motion (jWolf et al. 



19851). 



3.2. Power spectra of orbital elements 



The power spectrum, S{f), of an orbital element denoted by z(t) is defined as 



S{f) 



and is also useful to characterize the chaotic system. If one defines the autocorrelation 
function of z{t) as 

"^^M^ lim / z{t)z{t + T)dt, (9) 
its Fourier tr ansform is e qual to the power spectrum according to the Wiener- Khinchin 



theorem (e.g. 



Leichl 



19801 ): 



S{f) 



(10) 



Hence by investigating the power spectrum, we can acquire the knowledge of the temporal 
autocorrelation of the orbital element. 

The integrable systems have a discrete spectrum with peaks at the orbital periods 
provided appropriate variables are chosen. On the other hand, the chaotic systems have a 
rather featureless continuum spectrum. In particular, it is known that the nearly integrable 
systems show in general a power-law spectrum in their low frequency regime as 



s{f) oc r 



indicating a long-time correlation for the variables. 



In this paper, we discuss the power spectra of the orbital eccentricity, e^, for each planet 
in the system. We confirmed, however, that other orbital elements such as the semi-major 
axis give a similar result. 



- 11 - 



3.3. Distribution of induction periods 



It is known for the induction phenomena that the distribution of the induction periods, 
T, or the duration of quasi-regular motions generally obeys a power law, 



P(T) oc T-^ {13 > 0), 



(12) 



for large T (IBaouer fc BertschI Il990l ). This implies that there is no characteristic time 
scale for the induction period. This type of distribution can be obtained in the stagnant 
motion model by evoking a collection of tori, which have a self-simil ar distribution, in 



19891 ). 



the so-called stagnant layer around a KAM torus in the phase space (lAizawa et al. 
Fig. [1] illustrates schematically the phase space structure assumed in the stagnant motion 
model. The stagnant layer is filled with self-similarly distributed tori, which trap the 
system around them for a long time. Since the trapping time is scale free thanks to the 
self-similar distributions of the tori, the power law is obtained for the induction period, 
which is nothing but the trapping time. 

We will use this feature to see if the planetary motions of our concern are indeed the 
induction phenomenon. We expect that the duration of quasi-regular motions corresponds 
to the induction period. More precisely, we define the duration of the regular motion as the 
interval from the start of the integration of motion until the first close encounter between 
two planets. The close encounter means here that the pair of planets has the distance 
between them smaller than their Hill radius. It is known from the previous papers and 
confirmed in this paper that the transition fro m the regular motion t o the chaotic one 



19961 ). Incidentally, since 



occurs in general after the first close encounter (iChambers et al.l 
the ejection of one planet occurs rather soon after the encounter for the Jupiter system, we 
also study the distribution of the time from the encounter to the ejection. 
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4. Results 

4.1. Lyapunov exponents 

We first show in Fig. |2]the evolution of the orbital semi-major axis and eccentricity of 
each planet for model 11, which is representative of the proto-planet system. In this model, 
the first close encounter between two planets occurs at t = 16400yr. Before the encounter, 
the planetary motions are almost regular and the orbital elements remain unchanged 
essentially. After the encounter, on the other hand, they start to vary on short time scales. 
Then comes a period it ^ 23500 — 38500yr) when the semi-major axises do not change 
very much and the eccentricities for two planets vary rather monotonically in this period. 
Thereafter the orbital elements change in time violently again. Note, however, that the 
quasi-regular phases, although with much shorter periods, emerge intermittently. 

In Fig. |3l both the global and local Lyapunov exponents, Xgiobaiit) and \iocai{n,T), are 
displayed as a function of time for the same model as in Fig. [2l As mentioned earlier, the 
time interval r for the integration of the local Lyapunov exponent is set to be lOOyr, that is 
roughly 100 times the orbital period for the proto-planet system, whose innermost planet is 
initially located at lAU. It is clear that Xgiobai{t) monotonically decreases toward zero before 
the close encounter (A ~ 0.0005 just prior to the encounter), indicating that the motions 
are non-chaotic (or very weakly chaotic) in this phase. After the encounter, on the other 
hand, Xgiobaiit) increases drastically by about two orders of magnitude and fluctuates slowly 
around a constant value (~ 0.02) there after. For compari s on, w e also employ MEGNO, 



another indicator of chaos suggested by 



Cincotta &: Simd (120001 ) . to estimate the global 



Lyapunov exponent. The values derived from MEGNO agree with the ones obtained above 
within 1% typically. 

The local Lyapunov exponent also shows that the orbital motions are almost non- 
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chaotic before the encounter. In fact, the almost constant small value during this period is 
consistent with the evolution of the global Lyapunov exponent for the first lOOyr. It is also 
clear that the local Lyapunov exponent shows remarkable peaks rather intermittently after 
the close encounter. Moreover, it is found by the comparison between Figs. [2] and [3] that 
^locai takes a small constant value when the semi-major axis of each planet remains nearly 
constant in time. It is interesting to point out again that the eccentricities of planetary 
orbits are not zero and change in time during this period. This suggests that the planetary 
system is settled to a quasi-stable configuration what is different from the initial condition 
and has substantial orbital eccentricities. 

In Fig. m we show the orbital evolution of the Jupiter system (model 17). In this 
model, the first close encounter happens at t = 162416yr and one of the planets is ejected 
from the system at t = 239279yr. Just as in the proto-planet system shown in Fig. [21 
the semi-major axises remain almost unchanged in time and the eccentricities oscillate 
around the initial value with small amplitudes before the close encounter but they start to 
vary rapidly in time after the encounter. It is noted that the amplitudes of the variations 
are much larger for the Jupiter system than those for the proto-planet system, the fact 
responsible for the ejection of a planet in short times in the Jupiter system. 

In Fig. [5l the global and local Lyapunov exponents are shown for the same Jupiter 
system. Again we see the monotonic decrease of the exponent for the first ~ 20000yr. In 
this case, however, the exponent is then saturated and stays at a small but finite level. 
This reflects the fact that the system is close to integrable but still non-integrable. The 
global Lyapunov exponent increases quickly after the close encounter as in the proto-planet 
system. The local Lyapunov exponent obtained for every lOOyr, which is about 10 times 
the orbital period for the Jupiter system, shows some intermittent spikes after the close 
encounter although the interval from the close encounter to the ejection of a planet is rather 



14 



short. It is also confirmed that the violent variations of the orbital elements occur in the 
spikes of the local Lyapunov exponent and that when Xiocai has a small, nearly constant 
value, the semi-major axises are not changed very much whereas the eccentricities are 
non-zero and fluctuate rather slowly. The Lyapunov exponent increases toward the ejection 
of the planet in this case. 

The above-mentioned features in the orbital evolutions as well as in the Lyapunov 
exponents are common to all the models. (The behavior of the Lyapunov exponent close 
to the ejection of a planet is an exception and no clear trend can be seen). This suggests 
that the underlying structure in the phase space is not different very much from each other. 
In particular, the quasi-regular motions before the close encounter, which are very weakly 
chaotic at most, strongly suggest the existence of the KAM torus. This may be true even of 
the periods that occur intermittently after the close encounter, in which the local Lyapunov 
exponent returns to a small value and the planetary motions become quasi-regular again. 
Then the following picture is inferred: the phase space consists of a chaotic sea and KAM 
tori surrounded by a stagnant layer that consists of satellite tori. The phase space orbits 
go from one system of KAM torus and stagnant layer to another through the chaotic sea. 
When the phase space orbit is moving around one of the KAM tori, the local Lyapunov 
exponent takes a small constant value whereas it becomes spiky once the phase space orbit 
moves into the chaotic sea. 

This picture is also supported by the plot in Fig. |6]of the action (Jg) and angle (^o) 
variables of the non-perturbed Hamiltonian Hxep in Eq. dH]) • These variables are given as 



In the above equations, the angle variable ^o,i is called the mean anomaly and the so-called 




(13) 



6'o,i = Ui- Ci sinwj. 



(14) 
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eccentric anomaly, Ui, is defined as 




l + ei 




-tan—, 

i ^ 



(15) 



where Cj denotes the eccentricity and (pi the orbital phase of the i—ih. planet measured in the 
barycentric coordinates. In Fig. [6l we plot the action and angle variables for the innermost 
planet every ten steps. The behavior of the variables for other planets is essentially the 
same. 

One can recognize some horizontal bands and voids in the figure. These bands are 
regions, where the phase space orbit lingers, whereas the voids are passed through quickly. 
It is also seen that some of the bands undulate. The densest band pointed by an arrow 
at Jo.i = 1.98 X 10^° (gcm^/s)) in the figure corresponds to the KAM torus of the initial 
regular motion and its satellite tori, whereas other bands represent other KAM tori and 
their satellites, which are visited by the phase space orbit during the evolution. The 
bands with undulation correspond to the motion approaching the ejection of a planet. 
Since the Jo,i — ^o,i plane is filled by horizontal lines (or tori) everywhere uniformly if the 
perturbations are absent, the voids can be interpreted as chaotic regions produced by the 
perturbations to Hxep- 

Looking more closely, one finds that all the lines composing a band that corresponds 
to a quasi-regular motion are oscillating with finite amplitudes. These oscillations are 
studied by the Fourier analysis and the power spectrum for the initial quasi-regular motion 
is plotted in Fig. [71 It is found that the power spectrum obeys a power-law over a wide 
frequency range. This means that there is no characteristic frequency scale, the fact which 
seems to be consistent with the stagnant motion model, in which satellite tori are supposed 
to exist around a KAM torus with a fractal size distribution. 

In the following we will further look for evidence that the dynamics in these periods, 
that is, a relatively long quasi-regular motions followed by an abrupt transition to chaotic 
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motions can be interpreted as an induction phenomenon indeed. 

4.2. Power spectra of orbital eccentricities 

In Fig. IHl we show the power spectra of the orbital eccentricity of the innermost planet 
before and after the close encounter for model 11. This is a representative model for the 
proto-planet system and the one we employed to demonstrate the behavior of the Lyapunov 
exponents in Fig. [31 This model also has a merit that it has a relatively long duration up to 
the first close encounter. Although the data are quite noisy, they can be roughly fit by the 
power law, oc l//'^. It is apparent that the power-law indices are different between the two 
phases. We obtain v = 1.01 before the encounter whereas the spectrum becomes steeper 
with u = 1.85 thereafter. 

The same trend can be seen also for the Jupiter system, whose typical results are 
displayed in Fig. IHlfor model 17. Again the spectra are fit by the power law approximately 
both before and after the encounter and the spectral indices are u = 0.88 before the 
encounter and u = 1.85 thereafter. The spectrum given in the upper panel of Fig. [9] has 
finer features than that in the upper panel of Fig. [8] because pre-encounter phase is longer 
for the model in Fig. [8] than for that in Fig. [91 Note also that the post-encounter phase was 
computed up to the ejection of a planet in this case. 

In Tables [21 and [31 we summarize the power-law spectral indices of the arbitrarily 
chosen three realizations for each model both before the encounter (1st phase in the table) 
and after the encounter (2nd phase). Although we again employ the orbital eccentricity of 
the innermost planet, the orbits of outer planets behave similarly. Table [2 gives the results 
for the proto-planet system while Table [31 corresponds to the Jupiter system. For some of 
the models, the duration of the phase is too short to obtain the spectral index and "-" is 
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put instead of the spectral index for them. Ahhough the spectral indices vary substantially 
even among different realizations for the same model, it is clear that phase 2 has spectral 
indices clustered around v = 2 whereas phase 1 has smaller indices, < z/ < 1.69 , in 
general. The average spectral indices over all the models in Table |2] are 0.84 and 1.84 for 
the phases 1 and 2, respectively. The counter parts for the Jupiter system given in Table [3] 
are 0.91 before the encounter and 1.77 after the encounter. 

The power-law spectrum is a characteristic feature of the stagnant motions although 
the power-law indices are not specified by the theory. The time variations of the orbital 
elements are induced by the energy exchange among planets. The power-law oc 1//^ 
observed for the post-encounter phase indicates that the time evolution of the orbital 
eccentricity is a Brownian motion with the root mean square being oc t^^^. On the other 
hand, the pre-encounter phase has a smaller spectral index in general. Although the 
numbers of the planets in the system are different, our re sults on the growt h rate of 



Zhou et al 



(120071 ) ■ who studied 



eccentricity in the post-encounter phase are consistent with 
the system with 50 planets and claimed that the eccentricity of planetary orbits roughly 
evolves as (e^)^/^ oc t^/^ both before and after the encounter. As mentioned above, however, 
we found for our system with three planets different power-law indices before the encounter. 
The smaller power-law indices of e(t) in the pre-encounter phase in our models might 
correspond to the so-called fractional Brownian motion, for which the root mean square 
g2^/2 grows more slowly than for the ordinary Brownian motion ( iMandelbrot fc Van Ness 



19681 ). 
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4.3. Distribution of induction periods 

4.3.1. proto-planet system 

As mentioned in section I3.3[ the distributions of tlie durations of various pliases sliould 
be one of tlie key ingredients if tlie dynamics were to be interpreted as an induction 
plienomenon and described by tlie stagnant motion model. For the proto-planet systems 
we study the statistics of the duration, Tee, of the pre-encounter phase, where the orbits are 
nearly circular. We show in Fig. [10] the distribution of Tee for 10000 realizations of each 
model for the proto-planet systems listed in Table [U For the plots we employed 100 equal 
bins between the maximum and minimum values of Tee. For some models with small initial 
orbital separations, these two values are quite different and, as a result, there are some bins 
with very small populations near the longest Tee. 

Except for the model with the smallest A, the distribution has a peak, which shifts to 
longer times as A becomes larger. It is also clear from the log-log plot that the dispersion 
around the peak gets larger, too, as A becomes greater. We are particularly interested in 
the long time regime, where the stagnant motion model predicts power laws, which implies 
that there is no characteristic time scale for the trapping in the stagnant layer and is 
supposed to be a consequence of the self-similar distribution of tori in the stagnant layer. In 
the figure, the straight lines are the power-law fit to the long-time part of the distributions. 
We employ 30 data points down from the one with the longest Tee. The bins with Tee 
shorter than that at the peak or those with the population of less than 2% are discarded. 
If there are less than 30 data points that satisfy the criteria, all of them are used. The 
obtained power-law indices lie between —5.68 to —2.27 as shown in each panel, which are 
expected to reflect the difference in the phase space structures. 

It is obvious that the distributions are deviated from power laws both at the short 
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and long durations. The initial conditions prepared so that planets are initially in regular 
motions may lead to the underestimation of the duration time, since it might have cut 
short the earlier portion of the pre-encounter phase. Note, however, that there is no reason 
that we expect power laws for short time scales. As for the longer time scales, on the other 
hand, power laws are expected if the stagnant motion model can be applied. We suspect 
that the main reason for the deviation from the power law for the very long durations is 
that the number of realizations, that is 10000 for each model, are not large enough. In 
fact, only a small number of realizations are contributing to the longest time portion of the 
distributions, in which the deviation from the power law is remarkable. 

In Fig. [m we present a histogram in the A — Tee plane expressed in color for the 
number of cases in our model calculations. Connecting the peak, Tce,peaki for each A, we 
find the following relation: 

logT,e,peafc ^ 1.008A - 1.307, (16) 



which is very similar to what IChambers et al. 



(119961 ) found (a = 1.176,6 = -1.663 in 



their results) with a much smaller number (5) of realizations. It is clear from our results 
that the relation holds only for the durations corresponding to the peaks and, in fact, the 
durations for each A have a distribution as demonstrated above. Incidentally, the integrals 
of motion, the linear and angular momenta and total energy, are fixed i n pro ducing different 



realizations in this paper, which was not the case for 



Chambers et al 



f ll996h . 



4^.3.2. Jupiter system 

One of the characteristics for the dynamics of the Jupiter systems is that one of 
the planets are ejected from the system eventually. In addition to the durations of the 
pre-encounter phase, T^e-, we take also the statistics of the time from the first close encounter 
until the ejection of a planet, T^j — T^e, for the Jupiter systems. 
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Fig. [12] shows the distributions of Tee for models 17 (left panel) and 18 (right 
panel) for the Jupiter syst e m. N ote that model 18 is meant to mimic the models in 



Marzari fc Weidenschillingl ( 120021 ) and the initial condition is prepared differently from 
other models. The straight lines in the figure are the power-law fit to the long-duration part 
of the distributions and obtained just in the same manner as for the proto-planet systems 
(see Fig. [TU]) . It is seen again that the power- law distribution is a good approximation in 
this regime, which suggests that the quasi-regular motions before the close encounter in the 
Jupiter system can be also understood as the stagnant motion just as for the proto-planet 
system. The obtained power-law indices are —5.57 and —1.56 for models 17 and 18, 
respectively. 

Now we turn attention to the distributions of the time from the encounter to the 
ejection, T^j — Tee, which are given in Fig. [12] for models 17 (left panel) and 18 (right 
panel). It is clear from the figure that the long-time portions of both distributions are 
again approximated by the power laws with the indices of —2.65 and —2.67 for models 17 
and 18, respectively. This implies that there is no characteristic time scale for the duration, 
during which the phase space orbit of the system remains in the part of the phase space 
corresponding to the bound states of three planets. 

If the power-law distribution of Tee reflects the self-similar distribution of smaller tori 
in the stagnant layer of a KAM torus (see Fig. [1]) as claimed in the stagnant motion model, 
the power-law distribution of Tej — Tee might suggest a self-similar distribution of these 
KAM tori in the part of the phase space for the bound motions of three planets (see Fig 
for a schematic picture of the phase space). 
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5. Summary and Discussions 

In this paper we have numerically investigated the dynamics of the three-planet system 
and inferred its phase space structure from the obtained Lyapunov exponents, power 
spectra of orbital elements and distributions of induction periods based on the stagnant 
motion model. What we have found are the foUowings: 

1. The global and local Lyapunov exponents show that the system is almost non-chaotic 
until the first close encounter between two planets and it then turns into chaotic motions 
with intermittent non-chaotic periods. This suggests that the phase space consists of KAM 
tori surrounded by the stagnant layer and immersed in the chaotic sea. In fact, the dense 
bands are formed in the action-angle plot, corresponding to the quasi-regular motions. The 
lines composing a band are undulating with frequencies that obey a power-law and this 
may represent the motions around the satellite tori in the stagnant layer, which have a 
self-similar distribution. The phase space orbit goes from one system of KAM torus and 
stagnant layer to another through the chaotic sea. 

2. The power spectra of the orbital eccentricities of planets can be approximated 
by the power law, oc l/f^, in general. Such power-law spectra are known to be one of 
the characteristic features of the stagnant motions although the power-law index is not 
predicted by the theory. In our models the power-law index is ~ 1 for the pre-encounter 
phase whereas it becomes z/ ^ 2 after the encounter. The spectrum in the post-encounter 
phase is similar to that of the Brownian motions or the random walks. On the other hand, 
the spectrum in the pre-encounter phase might be originated from the fractional Brownian 
motions. 

3. The distributions of the duration of the pre-encounter phase that was referred to as 
the induction periods obey the power law in the long-duration part. The power- law indices 
are substantially different between models. It is stressed that the stagnant motion model 
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predicts the power law for the distribution of the induction periods as a consequence of 
the self-similar distribution of smaller tori in the stagnant layer around a KAM torus. The 
distributions are deviated from the power law in the short- duration part and has a peak in 
between. Connecting the peaks for various models with differ e nt ini tial orbital separations, 



we have obtained the relation similar to what 



Chambers et al. 



( 119961 ) found. It is also shown 



that the duration of the pre-encounter phase has actually a considerably broad distribution. 

4. For the Jupiter system, the distribution of the time from the first encounter to 
the ejection of a planet from the system also obeys a power law, which was not expected 
initially. From the analogy to the stagnant motion model, we might be able to infer the 
phase space structure as shown schematically in Fig. [TH Many KAM tori with its own 
stagnant layer and satellite tori in it are distributed self-similarly in the chaotic sea. 

Although we expect the phase space structure depicted in Fig. [H] is true both of the 
proto-planet system and the Jupiter system, the difference between them should be also 
mentioned. In general, the number of KAM tori in the phase space becomes smaller and 
the stagnant layers around them get thinner as the perturbation to the integrable system 
is greater. In the system of our concern, the perturbation is the interactions among the 
planets and hence it is larger for more massive planets. The pre-encounter phase is an 
exception, though. In this phase, the planets have nearly circular orbits separated by 
several Hill radii. Then the strength of the interactions between the planets depends only 
weakly on the planetary mass thanks to the definition of the Hill radius given in Eq. ([2]). 
These facts suggest that the KAM torus corresponding to the initial regular motion and its 
stagnant layer are robust and similar for the proto-planet and Jupiter systems while the 
number of other KAM tori is smaller and their stagnant layers are thinner for the Jupiter 
system than for the proto-planet system. This difference in the phase space structures is 
supposed to be responsible for the observed difference in the orbital evolutions of the two 
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systems: one of the planets is ejected in short times for the Jupiter system whereas no 
ejection occurs at least for ~ 10^ years in the proto-planet system. 

The results obtained in this paper appear to be consistent with our interpretation that 
the dynamics of three-planet system is a stagnant motion at least in the pre-encounter phase. 
It is also suggested from the results for the Jupiter system that even the post-encounter 
phase may be described by some extension of the stagnant motion model. It is true, 
however, that a more direct capture of satellite tori in the phase space is certainly desirable. 
Although we have attempted to do this with the so-called Poincare mapping, but in vain 
so far. We are afraid that the degree of freedom of our system is just too large to find 
an appropriate two dimensional section in the 12-dimensional phase space. Maybe other 
approaches should be pursued in the future work. In so doing, the number of realizations 
should be increased and other initial settings should be tried. Not to mention, we are also 
interested in how the results will change as the number of planets are varied. 
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Fig. 1. — A schematic illustration of the structure in the phase space proposed by the 
stagnant motion model. The KAM torus is surrounded by many smaller tori distributed in 
a self-similar manner in the stagnant layer. An exemplary evolution in the phase space is 
shown by a curve with an arrow. The system is trapped by the tori in the stagnant layer 
and stay there for a long time, showing quasi-regular motions, which are turned into chaotic 
motions when the phase-space orbit escapes out of the stagnant layer into the chaotic sea. 
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(lower panel) of each planet for model 11. 
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exponents for the same model as in Fig. [2J See the text for their definitions. 
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Fig. 8. — Power spectra S{f) of the orbital eccentricity of the innermost planet before (upper 
panel) and after (lower panel) the close encounter for model 11. The straight lines are the 
fit to the data. The power-law indices are 1.01 and 1.85 before and after the encounter, 
respectively. 
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Fig. 9. — Power spectra S{f) of the orbital eccentricity of the innermost planet before (upper 
panel) and after (lower panel) the close encounter for model 17. The straight lines are the 
fit to the data. The power-law indices are 0.88 and 1.85 before and after the encounter, 
respectively. 
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Fig. 10. — Normalized distributions of the duration of the pre-encounter phase, Tee, for 
models 1 to 16. The straight lines are the fit by power law to the data in the long-duration 
regime. 
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Fig. 11. — Histogram for models 1 through 16 in the plane of the initial orbital separation 
A and the duration of the pre-encounter phase, T^e- The number of cases is expressed in 
color. The blue dotted line is the fit to the pe ak (Eq. ( IT6l) ) and the green dashed line shows 



the A — relation in 



Chambers et al 



19961 ). 
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Fig. 12. — Normalized distributions of the duration of the pre-encounter phase, T^e, for 
models 17 (left panel) and 18 (right panel). The straight lines are the fit by power law to 
the data in the long-duration regime. 
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Fig. 13. — Normalized distributions of the time from the first close encounter up to the 
ejection of a planet, Tgj — Tee, for models 17 (left panel) and 18 (right panel). The straight 
lines are the fit by power law to the data in the long-time regime. 
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Fig. 14. — A schematic illustration of the phase space suggested by the power-law distribu- 
tions of the time from the encounter to the ejection of a planet for the Jupiter system. KAM 
tori of various sizes with a stagnant layer and smaller tori in it are distributed self-similarly 
in the chaotic sea. An exemplary phase space orbit is given by a curve with an arrow. 
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13 


10-' 


1 


-2.5902 


2 


6883 


-5.5261 


9.6417 


5 


3987 


(1 


8791 


X 


10-5) 


1 


2058 


X 


10-' 


(3 


7013 


X 


10- 


15 


14 


10-' 


1 


-2.5881 


2 


6894 


-5.7273 


9.9955 


5 


5988 


(1 


6665 


X 


10-5) 


1 


2052 


X 


10-' 


(3 


6969 


X 


10- 


15 


15 


10-' 


1 


-2.5860 


2 


6905 


-5.9282 


10.349 


5 


7988 


(1 


4844 


X 


10-5) 


1 


2047 


X 


10-' 


(3 


6916 


X 


10- 


15 


16 


10-' 


1 


-2.5840 


2 


6916 


-6.1290 


10.702 


5 


9989 


(1 


3274 


X 


10-5) 


1 


2040 


X 


10-' 


(3 


6874 


X 


10- 


15 


17 


10-^ 


5 


-4.0528 X 10^ 


6.9568 X 10* 


-2.7455 X 10^5 


5.2670 X 10^ 


3 


3968 


(3 


3889 


X 


10-*) 


1.0252 


X 


10-3 


(2.4813 


X 


10" 




18 


10-3 


5 


-3.9198 X 10^ 


7.1002 X lO* 


-2.4088 X 10^ 


5.8830 X 10^ 


^ < a2 


>= 7.2509 (8.5876 X 10-*) 


1.0138 


X 


10-3 


(2.4088 


X 


10" 





,<a3>= 9.4989 (1.1202 x 10-3) 



Note. — ai : the semi-major axis of the innermost planet. Etot '■ the total energy. L^tot '■ the z-component of the total angular momentum. Pxtot, Pytot '■ the 
X- and y-componcnts of the total linear momentum. < A >, tr^ : the mean and variance of the planetary orbital separations over 10000 realizations. < e >, o-g : 
the mean and variance of the planetary eccentricity over 10000 realizations. For each realization the mean eccentricity over three planets, e, is obtained and then 
< e > and a'^ are computed. 

< a2 >, T^^, < 03 >, fT^^ : the mean and variance of the semi- major axis of the second and third planets. 
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Table 2. The Power Spectral Index for the Proto-Planet Systems 



Model 


Orbital separation 




power- law index v 








sample 1 


sample 2 


sample 3 




A (Hill radius) 


1st phase 2nd phase 


1st phase 2nd phase 


1st phase 2nd phase 


1 


2.9970 


2.12 


1.85 


1.95 


2 


3.1973 


1.59 


2.06 


1.36 


3 


3.3975 


1.61 


2.38 


1.79 


4 


3.5977 


1.83 


1.59 


2.03 


5 


3.7979 


2.00 


1.70 


1.43 


6 


3.9980 


2.03 


1.52 


1.90 


7 


4.1982 


1.91 


2.00 


1.73 


8 


4.3983 


1.66 


2.01 


1.79 


9 


4.5984 


1.67 


1.66 


2.08 


10 


4.7985 


0.10 1.89 


0.65 2.01 


1.27 1.39 


11 


4.9986 


2.12 


1.01 1.85 


1.95 


12 


5.1986 


0.75 1.62 


1.92 


0.91 2.31 


13 


5.3987 


0.83 0.93 


0.93 1.90 


0.68 1.96 


14 


5.5988 


0.84 1.85 


1.25 1.97 


0.85 1.96 


15 


5.7988 


0.92 1.88 


0.92 1.93 


0.50 1.86 


16 


5.9989 


1.09 1.80 


0.65 1.84 


0.90 1.95 



Note. — The power spectra of the orbital eccentricity for the innermost planet are fit by the power law. 
Three realizations are arbitrarily chosen among 10000 realizations for each model and shown as samples 1, 2 
and 3. The pre- and post-encounter phases are treated separately and referred to as 1st and 2nd phases in 
the table, respectively. For some of the models, the duration of the phase is too short to obtain the power-law 
index and "-" is put instead of the index. The averages of v over all the models in this table are 0.84 before the 
encounter and 1.84 after the encounter. 
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Table 3. The Power Spectral Index for the Jupiter Systems 



Model 


Orbital separation 






power-law index u 










sample 1 


sample 2 


sample 3 






1st phase 


2nd phase 


1st phase 2nd phase 


1st phase 


2nd phase 


1 


< A >= 3.3968(Hill radius) 


0.88 


1.85 


1.69 1.55 


0.00 


1.72 


2 


< 02 >= 7.2509, < a:i >= 9.4989(AU) 


0.99 


1.87 


1.03 1.78 


0.99 


1.87 



Note. — The power spectra of the orbital eccentricity for the innermost planet are fit by the power law. The averages of u over the 
models in the table are 0.91 before the encounter and 1.77 after the encounter. 



